home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / src / haeberli / mpeg / jrevdct.c < prev    next >
C/C++ Source or Header  |  1994-08-01  |  44KB  |  1,462 lines
  1. /*
  2.  * jrevdct.c
  3.  *
  4.  * Copyright (C) 1991, 1992, Thomas G. Lane.
  5.  * This file is part of the Independent JPEG Group's software.
  6.  * For conditions of distribution and use, see the accompanying README file.
  7.  *
  8.  * This file contains the basic inverse-DCT transformation subroutine.
  9.  *
  10.  * This implementation is based on an algorithm described in
  11.  *   C. Loeffler, A. Ligtenberg and G. Moschytz, "Practical Fast 1-D DCT
  12.  *   Algorithms with 11 Multiplications", Proc. Int'l. Conf. on Acoustics,
  13.  *   Speech, and Signal Processing 1989 (ICASSP '89), pp. 988-991.
  14.  * The primary algorithm described there uses 11 multiplies and 29 adds.
  15.  * We use their alternate method with 12 multiplies and 32 adds.
  16.  * The advantage of this method is that no data path contains more than one
  17.  * multiplication; this allows a very simple and accurate implementation in
  18.  * scaled fixed-point arithmetic, with a minimal number of shifts.
  19.  * 
  20.  * I've made lots of modifications to attempt to take advantage of the
  21.  * sparse nature of the DCT matrices we're getting.  Although the logic
  22.  * is cumbersome, it's straightforward and the resulting code is much
  23.  * faster.
  24.  *
  25.  * A better way to do this would be to pass in the DCT block as a sparse
  26.  * matrix, perhaps with the difference cases encoded.
  27.  */
  28.  
  29. #include <string.h>
  30. #include "video.h"
  31. #include "proto.h"
  32.  
  33. #define GLOBAL            /* a function referenced thru EXTERNs */
  34.   
  35. /* We assume that right shift corresponds to signed division by 2 with
  36.  * rounding towards minus infinity.  This is correct for typical "arithmetic
  37.  * shift" instructions that shift in copies of the sign bit.  But some
  38.  * C compilers implement >> with an unsigned shift.  For these machines you
  39.  * must define RIGHT_SHIFT_IS_UNSIGNED.
  40.  * RIGHT_SHIFT provides a proper signed right shift of an INT32 quantity.
  41.  * It is only applied with constant shift counts.  SHIFT_TEMPS must be
  42.  * included in the variables of any routine using RIGHT_SHIFT.
  43.  */
  44.   
  45. #ifdef RIGHT_SHIFT_IS_UNSIGNED
  46. #define SHIFT_TEMPS    INT32 shift_temp;
  47. #define RIGHT_SHIFT(x,shft)  \
  48.     ((shift_temp = (x)) < 0 ? \
  49.      (shift_temp >> (shft)) | ((~((INT32) 0)) << (32-(shft))) : \
  50.      (shift_temp >> (shft)))
  51. #else
  52. #define SHIFT_TEMPS
  53. #define RIGHT_SHIFT(x,shft)    ((x) >> (shft))
  54. #endif
  55.  
  56. /*
  57.  * This routine is specialized to the case DCTSIZE = 8.
  58.  */
  59.  
  60. #if DCTSIZE != 8
  61.   Sorry, this code only copes with 8x8 DCTs. /* deliberate syntax err */
  62. #endif
  63.  
  64.  
  65. /*
  66.  * A 2-D IDCT can be done by 1-D IDCT on each row followed by 1-D IDCT
  67.  * on each column.  Direct algorithms are also available, but they are
  68.  * much more complex and seem not to be any faster when reduced to code.
  69.  *
  70.  * The poop on this scaling stuff is as follows:
  71.  *
  72.  * Each 1-D IDCT step produces outputs which are a factor of sqrt(N)
  73.  * larger than the true IDCT outputs.  The final outputs are therefore
  74.  * a factor of N larger than desired; since N=8 this can be cured by
  75.  * a simple right shift at the end of the algorithm.  The advantage of
  76.  * this arrangement is that we save two multiplications per 1-D IDCT,
  77.  * because the y0 and y4 inputs need not be divided by sqrt(N).
  78.  *
  79.  * We have to do addition and subtraction of the integer inputs, which
  80.  * is no problem, and multiplication by fractional constants, which is
  81.  * a problem to do in integer arithmetic.  We multiply all the constants
  82.  * by CONST_SCALE and convert them to integer constants (thus retaining
  83.  * CONST_BITS bits of precision in the constants).  After doing a
  84.  * multiplication we have to divide the product by CONST_SCALE, with proper
  85.  * rounding, to produce the correct output.  This division can be done
  86.  * cheaply as a right shift of CONST_BITS bits.  We postpone shifting
  87.  * as long as possible so that partial sums can be added together with
  88.  * full fractional precision.
  89.  *
  90.  * The outputs of the first pass are scaled up by PASS1_BITS bits so that
  91.  * they are represented to better-than-integral precision.  These outputs
  92.  * require BITS_IN_JSAMPLE + PASS1_BITS + 3 bits; this fits in a 16-bit word
  93.  * with the recommended scaling.  (To scale up 12-bit sample data further, an
  94.  * intermediate INT32 array would be needed.)
  95.  *
  96.  * To avoid overflow of the 32-bit intermediate results in pass 2, we must
  97.  * have BITS_IN_JSAMPLE + CONST_BITS + PASS1_BITS <= 26.  Error analysis
  98.  * shows that the values given below are the most effective.
  99.  */
  100.  
  101. #ifdef EIGHT_BIT_SAMPLES
  102. #define PASS1_BITS  2
  103. #else
  104. #define PASS1_BITS  1        /* lose a little precision to avoid overflow */
  105. #endif
  106.  
  107. #define ONE    ((INT32) 1)
  108.  
  109. #define CONST_SCALE (ONE << CONST_BITS)
  110.  
  111. /* Convert a positive real constant to an integer scaled by CONST_SCALE.
  112.  * IMPORTANT: if your compiler doesn't do this arithmetic at compile time,
  113.  * you will pay a significant penalty in run time.  In that case, figure
  114.  * the correct integer constant values and insert them by hand.
  115.  */
  116.  
  117. #define FIX(x)    ((INT32) ((x) * CONST_SCALE + 0.5))
  118.  
  119. /* Descale and correctly round an INT32 value that's scaled by N bits.
  120.  * We assume RIGHT_SHIFT rounds towards minus infinity, so adding
  121.  * the fudge factor is correct for either sign of X.
  122.  */
  123.  
  124. #define DESCALE(x,n)  RIGHT_SHIFT((x) + (ONE << ((n)-1)), n)
  125.  
  126. /* Multiply an INT32 variable by an INT32 constant to yield an INT32 result.
  127.  * For 8-bit samples with the recommended scaling, all the variable
  128.  * and constant values involved are no more than 16 bits wide, so a
  129.  * 16x16->32 bit multiply can be used instead of a full 32x32 multiply;
  130.  * this provides a useful speedup on many machines.
  131.  * There is no way to specify a 16x16->32 multiply in portable C, but
  132.  * some C compilers will do the right thing if you provide the correct
  133.  * combination of casts.
  134.  * NB: for 12-bit samples, a full 32-bit multiplication will be needed.
  135.  */
  136.  
  137. #ifdef EIGHT_BIT_SAMPLES
  138. #ifdef SHORTxSHORT_32        /* may work if 'int' is 32 bits */
  139. #define MULTIPLY(var,const)  (((INT16) (var)) * ((INT16) (const)))
  140. #endif
  141. #ifdef SHORTxLCONST_32        /* known to work with Microsoft C 6.0 */
  142. #define MULTIPLY(var,const)  (((INT16) (var)) * ((INT32) (const)))
  143. #endif
  144. #endif
  145.  
  146. #ifndef MULTIPLY        /* default definition */
  147. #define MULTIPLY(var,const)  ((var) * (const))
  148. #endif
  149.  
  150. /* Precomputed idct value arrays. */
  151.  
  152. static DCTELEM PreIDCT[64][64];
  153.  
  154. /* Pre compute singleton coefficient IDCT values. */
  155. void
  156. init_pre_idct() {
  157.   int i;
  158.   void j_rev_dct();
  159.  
  160.   for (i=0; i<64; i++) {
  161.     memset((char *) PreIDCT[i], 0, 64*sizeof(DCTELEM));
  162.     PreIDCT[i][i] = 2048;
  163.     j_rev_dct(PreIDCT[i]);
  164.   }
  165. }
  166.  
  167. #ifndef ORIG_DCT
  168.   
  169.  
  170. /*
  171.  * Perform the inverse DCT on one block of coefficients.
  172.  */
  173.  
  174. void
  175. j_rev_dct_sparse (data, pos)
  176.      DCTBLOCK data;
  177.      int pos;
  178. {
  179.   register DCTELEM *dataptr;
  180.   short int val;
  181.   DCTELEM *ndataptr;
  182.   int scale, coeff, rr;
  183.   register int *dp;
  184.   register int v;
  185.  
  186.   /* If DC Coefficient. */
  187.   
  188.   if (pos == 0) {
  189.     dp = (int *)data;
  190.     v = *data;
  191.     /* Compute 32 bit value to assign.  This speeds things up a bit */
  192.     if (v < 0) val = (v-3)>>3;
  193.     else val = (v+4)>>3;
  194.     v = val | (val << 16);
  195.     dp[0] = v;      dp[1] = v;      dp[2] = v;      dp[3] = v;
  196.     dp[4] = v;      dp[5] = v;      dp[6] = v;      dp[7] = v;
  197.     dp[8] = v;      dp[9] = v;      dp[10] = v;      dp[11] = v;
  198.     dp[12] = v;      dp[13] = v;      dp[14] = v;      dp[15] = v;
  199.     dp[16] = v;      dp[17] = v;      dp[18] = v;      dp[19] = v;
  200.     dp[20] = v;      dp[21] = v;      dp[22] = v;      dp[23] = v;
  201.     dp[24] = v;      dp[25] = v;      dp[26] = v;      dp[27] = v;
  202.     dp[28] = v;      dp[29] = v;      dp[30] = v;      dp[31] = v;
  203.     return;
  204.   }
  205.   
  206.   /* Some other coefficient. */
  207.   dataptr = (DCTELEM *)data;
  208.   coeff = dataptr[pos];
  209.   ndataptr = PreIDCT[pos];
  210.  
  211.   for (rr=0; rr<4; rr++) {
  212.     dataptr[0] = (ndataptr[0] * coeff) >> (CONST_BITS-2);
  213.     dataptr[1] = (ndataptr[1] * coeff) >> (CONST_BITS-2);
  214.     dataptr[2] = (ndataptr[2] * coeff) >> (CONST_BITS-2);
  215.     dataptr[3] = (ndataptr[3] * coeff) >> (CONST_BITS-2);
  216.     dataptr[4] = (ndataptr[4] * coeff) >> (CONST_BITS-2);
  217.     dataptr[5] = (ndataptr[5] * coeff) >> (CONST_BITS-2);
  218.     dataptr[6] = (ndataptr[6] * coeff) >> (CONST_BITS-2);
  219.     dataptr[7] = (ndataptr[7] * coeff) >> (CONST_BITS-2);
  220.     dataptr[8] = (ndataptr[8] * coeff) >> (CONST_BITS-2);
  221.     dataptr[9] = (ndataptr[9] * coeff) >> (CONST_BITS-2);
  222.     dataptr[10] = (ndataptr[10] * coeff) >> (CONST_BITS-2);
  223.     dataptr[11] = (ndataptr[11] * coeff) >> (CONST_BITS-2);
  224.     dataptr[12] = (ndataptr[12] * coeff) >> (CONST_BITS-2);
  225.     dataptr[13] = (ndataptr[13] * coeff) >> (CONST_BITS-2);
  226.     dataptr[14] = (ndataptr[14] * coeff) >> (CONST_BITS-2);
  227.     dataptr[15] = (ndataptr[15] * coeff) >> (CONST_BITS-2);
  228.     dataptr += 16;
  229.     ndataptr += 16;
  230.   }
  231.   return;
  232. }
  233.  
  234.  
  235. void
  236. j_rev_dct (data)
  237.      DCTBLOCK data;
  238. {
  239.   INT32 tmp0, tmp1, tmp2, tmp3;
  240.   INT32 tmp10, tmp11, tmp12, tmp13;
  241.   INT32 z1, z2, z3, z4, z5;
  242.   INT32 d0, d1, d2, d3, d4, d5, d6, d7;
  243.   register DCTELEM *dataptr;
  244.   int rowctr;
  245.   SHIFT_TEMPS
  246.    
  247.   /* Pass 1: process rows. */
  248.   /* Note results are scaled up by sqrt(8) compared to a true IDCT; */
  249.   /* furthermore, we scale the results by 2**PASS1_BITS. */
  250.  
  251.   dataptr = data;
  252.  
  253.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  254.     /* Due to quantization, we will usually find that many of the input
  255.      * coefficients are zero, especially the AC terms.  We can exploit this
  256.      * by short-circuiting the IDCT calculation for any row in which all
  257.      * the AC terms are zero.  In that case each output is equal to the
  258.      * DC coefficient (with scale factor as needed).
  259.      * With typical images and quantization tables, half or more of the
  260.      * row DCT calculations can be simplified this way.
  261.      */
  262.  
  263.     register int *idataptr = (int*)dataptr;
  264.     d0 = dataptr[0];
  265.     d1 = dataptr[1];
  266.     if ((d1 == 0) && (idataptr[1] | idataptr[2] | idataptr[3]) == 0) {
  267.       /* AC terms all zero */
  268.       if (d0) {
  269.       /* Compute a 32 bit value to assign. */
  270.       DCTELEM dcval = (DCTELEM) (d0 << PASS1_BITS);
  271.       register int v = (dcval & 0xffff) | ((dcval << 16) & 0xffff0000);
  272.       
  273.       idataptr[0] = v;
  274.       idataptr[1] = v;
  275.       idataptr[2] = v;
  276.       idataptr[3] = v;
  277.       }
  278.       
  279.       dataptr += DCTSIZE;    /* advance pointer to next row */
  280.       continue;
  281.     }
  282.     d2 = dataptr[2];
  283.     d3 = dataptr[3];
  284.     d4 = dataptr[4];
  285.     d5 = dataptr[5];
  286.     d6 = dataptr[6];
  287.     d7 = dataptr[7];
  288.  
  289.     /* Even part: reverse the even part of the forward DCT. */
  290.     /* The rotator is sqrt(2)*c(-6). */
  291.     if (d6) {
  292.     if (d4) {
  293.         if (d2) {
  294.         if (d0) {
  295.             /* d0 != 0, d2 != 0, d4 != 0, d6 != 0 */
  296.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  297.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  298.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  299.  
  300.             tmp0 = (d0 + d4) << CONST_BITS;
  301.             tmp1 = (d0 - d4) << CONST_BITS;
  302.  
  303.             tmp10 = tmp0 + tmp3;
  304.             tmp13 = tmp0 - tmp3;
  305.             tmp11 = tmp1 + tmp2;
  306.             tmp12 = tmp1 - tmp2;
  307.         } else {
  308.             /* d0 == 0, d2 != 0, d4 != 0, d6 != 0 */
  309.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  310.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  311.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  312.  
  313.             tmp0 = d4 << CONST_BITS;
  314.  
  315.             tmp10 = tmp0 + tmp3;
  316.             tmp13 = tmp0 - tmp3;
  317.             tmp11 = tmp2 - tmp0;
  318.             tmp12 = -(tmp0 + tmp2);
  319.         }
  320.         } else {
  321.         if (d0) {
  322.             /* d0 != 0, d2 == 0, d4 != 0, d6 != 0 */
  323.             tmp2 = MULTIPLY(d6, - FIX(1.306562965));
  324.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  325.  
  326.             tmp0 = (d0 + d4) << CONST_BITS;
  327.             tmp1 = (d0 - d4) << CONST_BITS;
  328.  
  329.             tmp10 = tmp0 + tmp3;
  330.             tmp13 = tmp0 - tmp3;
  331.             tmp11 = tmp1 + tmp2;
  332.             tmp12 = tmp1 - tmp2;
  333.         } else {
  334.             /* d0 == 0, d2 == 0, d4 != 0, d6 != 0 */
  335.             tmp2 = MULTIPLY(d6, -FIX(1.306562965));
  336.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  337.  
  338.             tmp0 = d4 << CONST_BITS;
  339.  
  340.             tmp10 = tmp0 + tmp3;
  341.             tmp13 = tmp0 - tmp3;
  342.             tmp11 = tmp2 - tmp0;
  343.             tmp12 = -(tmp0 + tmp2);
  344.         }
  345.         }
  346.     } else {
  347.         if (d2) {
  348.         if (d0) {
  349.             /* d0 != 0, d2 != 0, d4 == 0, d6 != 0 */
  350.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  351.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  352.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  353.  
  354.             tmp0 = d0 << CONST_BITS;
  355.  
  356.             tmp10 = tmp0 + tmp3;
  357.             tmp13 = tmp0 - tmp3;
  358.             tmp11 = tmp0 + tmp2;
  359.             tmp12 = tmp0 - tmp2;
  360.         } else {
  361.             /* d0 == 0, d2 != 0, d4 == 0, d6 != 0 */
  362.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  363.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  364.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  365.  
  366.             tmp10 = tmp3;
  367.             tmp13 = -tmp3;
  368.             tmp11 = tmp2;
  369.             tmp12 = -tmp2;
  370.         }
  371.         } else {
  372.         if (d0) {
  373.             /* d0 != 0, d2 == 0, d4 == 0, d6 != 0 */
  374.             tmp2 = MULTIPLY(d6, - FIX(1.306562965));
  375.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  376.  
  377.             tmp0 = d0 << CONST_BITS;
  378.  
  379.             tmp10 = tmp0 + tmp3;
  380.             tmp13 = tmp0 - tmp3;
  381.             tmp11 = tmp0 + tmp2;
  382.             tmp12 = tmp0 - tmp2;
  383.         } else {
  384.             /* d0 == 0, d2 == 0, d4 == 0, d6 != 0 */
  385.             tmp2 = MULTIPLY(d6, - FIX(1.306562965));
  386.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  387.  
  388.             tmp10 = tmp3;
  389.             tmp13 = -tmp3;
  390.             tmp11 = tmp2;
  391.             tmp12 = -tmp2;
  392.         }
  393.         }
  394.     }
  395.     } else {
  396.     if (d4) {
  397.         if (d2) {
  398.         if (d0) {
  399.             /* d0 != 0, d2 != 0, d4 != 0, d6 == 0 */
  400.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  401.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  402.  
  403.             tmp0 = (d0 + d4) << CONST_BITS;
  404.             tmp1 = (d0 - d4) << CONST_BITS;
  405.  
  406.             tmp10 = tmp0 + tmp3;
  407.             tmp13 = tmp0 - tmp3;
  408.             tmp11 = tmp1 + tmp2;
  409.             tmp12 = tmp1 - tmp2;
  410.         } else {
  411.             /* d0 == 0, d2 != 0, d4 != 0, d6 == 0 */
  412.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  413.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  414.  
  415.             tmp0 = d4 << CONST_BITS;
  416.  
  417.             tmp10 = tmp0 + tmp3;
  418.             tmp13 = tmp0 - tmp3;
  419.             tmp11 = tmp2 - tmp0;
  420.             tmp12 = -(tmp0 + tmp2);
  421.         }
  422.         } else {
  423.         if (d0) {
  424.             /* d0 != 0, d2 == 0, d4 != 0, d6 == 0 */
  425.             tmp10 = tmp13 = (d0 + d4) << CONST_BITS;
  426.             tmp11 = tmp12 = (d0 - d4) << CONST_BITS;
  427.         } else {
  428.             /* d0 == 0, d2 == 0, d4 != 0, d6 == 0 */
  429.             tmp10 = tmp13 = d4 << CONST_BITS;
  430.             tmp11 = tmp12 = -tmp10;
  431.         }
  432.         }
  433.     } else {
  434.         if (d2) {
  435.         if (d0) {
  436.             /* d0 != 0, d2 != 0, d4 == 0, d6 == 0 */
  437.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  438.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  439.  
  440.             tmp0 = d0 << CONST_BITS;
  441.  
  442.             tmp10 = tmp0 + tmp3;
  443.             tmp13 = tmp0 - tmp3;
  444.             tmp11 = tmp0 + tmp2;
  445.             tmp12 = tmp0 - tmp2;
  446.         } else {
  447.             /* d0 == 0, d2 != 0, d4 == 0, d6 == 0 */
  448.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  449.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  450.  
  451.             tmp10 = tmp3;
  452.             tmp13 = -tmp3;
  453.             tmp11 = tmp2;
  454.             tmp12 = -tmp2;
  455.         }
  456.         } else {
  457.         if (d0) {
  458.             /* d0 != 0, d2 == 0, d4 == 0, d6 == 0 */
  459.             tmp10 = tmp13 = tmp11 = tmp12 = d0 << CONST_BITS;
  460.         } else {
  461.             /* d0 == 0, d2 == 0, d4 == 0, d6 == 0 */
  462.             tmp10 = tmp13 = tmp11 = tmp12 = 0;
  463.         }
  464.         }
  465.     }
  466.     }
  467.  
  468.  
  469.     /* Odd part per figure 8; the matrix is unitary and hence its
  470.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  471.      */
  472.  
  473.     if (d7) {
  474.     if (d5) {
  475.         if (d3) {
  476.         if (d1) {
  477.             /* d1 != 0, d3 != 0, d5 != 0, d7 != 0 */
  478.             z1 = d7 + d1;
  479.             z2 = d5 + d3;
  480.             z3 = d7 + d3;
  481.             z4 = d5 + d1;
  482.             z5 = MULTIPLY(z3 + z4, FIX(1.175875602));
  483.             
  484.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  485.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  486.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  487.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  488.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  489.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  490.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  491.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  492.             
  493.             z3 += z5;
  494.             z4 += z5;
  495.             
  496.             tmp0 += z1 + z3;
  497.             tmp1 += z2 + z4;
  498.             tmp2 += z2 + z3;
  499.             tmp3 += z1 + z4;
  500.         } else {
  501.             /* d1 == 0, d3 != 0, d5 != 0, d7 != 0 */
  502.             z1 = d7;
  503.             z2 = d5 + d3;
  504.             z3 = d7 + d3;
  505.             z5 = MULTIPLY(z3 + d5, FIX(1.175875602));
  506.             
  507.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  508.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  509.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  510.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  511.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  512.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  513.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  514.             
  515.             z3 += z5;
  516.             z4 += z5;
  517.             
  518.             tmp0 += z1 + z3;
  519.             tmp1 += z2 + z4;
  520.             tmp2 += z2 + z3;
  521.             tmp3 = z1 + z4;
  522.         }
  523.         } else {
  524.         if (d1) {
  525.             /* d1 != 0, d3 == 0, d5 != 0, d7 != 0 */
  526.             z1 = d7 + d1;
  527.             z2 = d5;
  528.             z3 = d7;
  529.             z4 = d5 + d1;
  530.             z5 = MULTIPLY(z3 + z4, FIX(1.175875602));
  531.             
  532.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  533.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  534.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  535.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  536.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  537.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  538.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  539.             
  540.             z3 += z5;
  541.             z4 += z5;
  542.             
  543.             tmp0 += z1 + z3;
  544.             tmp1 += z2 + z4;
  545.             tmp2 = z2 + z3;
  546.             tmp3 += z1 + z4;
  547.         } else {
  548.             /* d1 == 0, d3 == 0, d5 != 0, d7 != 0 */
  549.             tmp0 = MULTIPLY(d7, - FIX(0.601344887)); 
  550.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  551.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  552.             tmp1 = MULTIPLY(d5, - FIX(0.509795578));
  553.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  554.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  555.             z5 = MULTIPLY(d5 + d7, FIX(1.175875602));
  556.             
  557.             z3 += z5;
  558.             z4 += z5;
  559.             
  560.             tmp0 += z3;
  561.             tmp1 += z4;
  562.             tmp2 = z2 + z3;
  563.             tmp3 = z1 + z4;
  564.         }
  565.         }
  566.     } else {
  567.         if (d3) {
  568.         if (d1) {
  569.             /* d1 != 0, d3 != 0, d5 == 0, d7 != 0 */
  570.             z1 = d7 + d1;
  571.             z3 = d7 + d3;
  572.             z5 = MULTIPLY(z3 + d1, FIX(1.175875602));
  573.             
  574.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  575.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  576.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  577.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  578.             z2 = MULTIPLY(d3, - FIX(2.562915447));
  579.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  580.             z4 = MULTIPLY(d1, - FIX(0.390180644));
  581.             
  582.             z3 += z5;
  583.             z4 += z5;
  584.             
  585.             tmp0 += z1 + z3;
  586.             tmp1 = z2 + z4;
  587.             tmp2 += z2 + z3;
  588.             tmp3 += z1 + z4;
  589.         } else {
  590.             /* d1 == 0, d3 != 0, d5 == 0, d7 != 0 */
  591.             z3 = d7 + d3;
  592.             
  593.             tmp0 = MULTIPLY(d7, - FIX(0.601344887)); 
  594.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  595.             tmp2 = MULTIPLY(d3, FIX(0.509795579));
  596.             z2 = MULTIPLY(d3, - FIX(2.562915447));
  597.             z5 = MULTIPLY(z3, FIX(1.175875602));
  598.             z3 = MULTIPLY(z3, - FIX(0.785694958));
  599.             
  600.             tmp0 += z3;
  601.             tmp1 = z2 + z5;
  602.             tmp2 += z3;
  603.             tmp3 = z1 + z5;
  604.         }
  605.         } else {
  606.         if (d1) {
  607.             /* d1 != 0, d3 == 0, d5 == 0, d7 != 0 */
  608.             z1 = d7 + d1;
  609.             z5 = MULTIPLY(z1, FIX(1.175875602));
  610.  
  611.             z1 = MULTIPLY(z1, FIX(0.275899379));
  612.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  613.             tmp0 = MULTIPLY(d7, - FIX(1.662939224)); 
  614.             z4 = MULTIPLY(d1, - FIX(0.390180644));
  615.             tmp3 = MULTIPLY(d1, FIX(1.111140466));
  616.  
  617.             tmp0 += z1;
  618.             tmp1 = z4 + z5;
  619.             tmp2 = z3 + z5;
  620.             tmp3 += z1;
  621.         } else {
  622.             /* d1 == 0, d3 == 0, d5 == 0, d7 != 0 */
  623.             tmp0 = MULTIPLY(d7, - FIX(1.387039845));
  624.             tmp1 = MULTIPLY(d7, FIX(1.175875602));
  625.             tmp2 = MULTIPLY(d7, - FIX(0.785694958));
  626.             tmp3 = MULTIPLY(d7, FIX(0.275899379));
  627.         }
  628.         }
  629.     }
  630.     } else {
  631.     if (d5) {
  632.         if (d3) {
  633.         if (d1) {
  634.             /* d1 != 0, d3 != 0, d5 != 0, d7 == 0 */
  635.             z2 = d5 + d3;
  636.             z4 = d5 + d1;
  637.             z5 = MULTIPLY(d3 + z4, FIX(1.175875602));
  638.             
  639.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  640.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  641.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  642.             z1 = MULTIPLY(d1, - FIX(0.899976223));
  643.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  644.             z3 = MULTIPLY(d3, - FIX(1.961570560));
  645.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  646.             
  647.             z3 += z5;
  648.             z4 += z5;
  649.             
  650.             tmp0 = z1 + z3;
  651.             tmp1 += z2 + z4;
  652.             tmp2 += z2 + z3;
  653.             tmp3 += z1 + z4;
  654.         } else {
  655.             /* d1 == 0, d3 != 0, d5 != 0, d7 == 0 */
  656.             z2 = d5 + d3;
  657.             
  658.             z5 = MULTIPLY(z2, FIX(1.175875602));
  659.             tmp1 = MULTIPLY(d5, FIX(1.662939225));
  660.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  661.             z2 = MULTIPLY(z2, - FIX(1.387039845));
  662.             tmp2 = MULTIPLY(d3, FIX(1.111140466));
  663.             z3 = MULTIPLY(d3, - FIX(1.961570560));
  664.             
  665.             tmp0 = z3 + z5;
  666.             tmp1 += z2;
  667.             tmp2 += z2;
  668.             tmp3 = z4 + z5;
  669.         }
  670.         } else {
  671.         if (d1) {
  672.             /* d1 != 0, d3 == 0, d5 != 0, d7 == 0 */
  673.             z4 = d5 + d1;
  674.             
  675.             z5 = MULTIPLY(z4, FIX(1.175875602));
  676.             z1 = MULTIPLY(d1, - FIX(0.899976223));
  677.             tmp3 = MULTIPLY(d1, FIX(0.601344887));
  678.             tmp1 = MULTIPLY(d5, - FIX(0.509795578));
  679.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  680.             z4 = MULTIPLY(z4, FIX(0.785694958));
  681.             
  682.             tmp0 = z1 + z5;
  683.             tmp1 += z4;
  684.             tmp2 = z2 + z5;
  685.             tmp3 += z4;
  686.         } else {
  687.             /* d1 == 0, d3 == 0, d5 != 0, d7 == 0 */
  688.             tmp0 = MULTIPLY(d5, FIX(1.175875602));
  689.             tmp1 = MULTIPLY(d5, FIX(0.275899380));
  690.             tmp2 = MULTIPLY(d5, - FIX(1.387039845));
  691.             tmp3 = MULTIPLY(d5, FIX(0.785694958));
  692.         }
  693.         }
  694.     } else {
  695.         if (d3) {
  696.         if (d1) {
  697.             /* d1 != 0, d3 != 0, d5 == 0, d7 == 0 */
  698.             z5 = d1 + d3;
  699.             tmp3 = MULTIPLY(d1, FIX(0.211164243));
  700.             tmp2 = MULTIPLY(d3, - FIX(1.451774981));
  701.             z1 = MULTIPLY(d1, FIX(1.061594337));
  702.             z2 = MULTIPLY(d3, - FIX(2.172734803));
  703.             z4 = MULTIPLY(z5, FIX(0.785694958));
  704.             z5 = MULTIPLY(z5, FIX(1.175875602));
  705.             
  706.             tmp0 = z1 - z4;
  707.             tmp1 = z2 + z4;
  708.             tmp2 += z5;
  709.             tmp3 += z5;
  710.         } else {
  711.             /* d1 == 0, d3 != 0, d5 == 0, d7 == 0 */
  712.             tmp0 = MULTIPLY(d3, - FIX(0.785694958));
  713.             tmp1 = MULTIPLY(d3, - FIX(1.387039845));
  714.             tmp2 = MULTIPLY(d3, - FIX(0.275899379));
  715.             tmp3 = MULTIPLY(d3, FIX(1.175875602));
  716.         }
  717.         } else {
  718.         if (d1) {
  719.             /* d1 != 0, d3 == 0, d5 == 0, d7 == 0 */
  720.             tmp0 = MULTIPLY(d1, FIX(0.275899379));
  721.             tmp1 = MULTIPLY(d1, FIX(0.785694958));
  722.             tmp2 = MULTIPLY(d1, FIX(1.175875602));
  723.             tmp3 = MULTIPLY(d1, FIX(1.387039845));
  724.         } else {
  725.             /* d1 == 0, d3 == 0, d5 == 0, d7 == 0 */
  726.             tmp0 = tmp1 = tmp2 = tmp3 = 0;
  727.         }
  728.         }
  729.     }
  730.     }
  731.  
  732.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  733.  
  734.     dataptr[0] = (DCTELEM) DESCALE(tmp10 + tmp3, CONST_BITS-PASS1_BITS);
  735.     dataptr[7] = (DCTELEM) DESCALE(tmp10 - tmp3, CONST_BITS-PASS1_BITS);
  736.     dataptr[1] = (DCTELEM) DESCALE(tmp11 + tmp2, CONST_BITS-PASS1_BITS);
  737.     dataptr[6] = (DCTELEM) DESCALE(tmp11 - tmp2, CONST_BITS-PASS1_BITS);
  738.     dataptr[2] = (DCTELEM) DESCALE(tmp12 + tmp1, CONST_BITS-PASS1_BITS);
  739.     dataptr[5] = (DCTELEM) DESCALE(tmp12 - tmp1, CONST_BITS-PASS1_BITS);
  740.     dataptr[3] = (DCTELEM) DESCALE(tmp13 + tmp0, CONST_BITS-PASS1_BITS);
  741.     dataptr[4] = (DCTELEM) DESCALE(tmp13 - tmp0, CONST_BITS-PASS1_BITS);
  742.  
  743.     dataptr += DCTSIZE;        /* advance pointer to next row */
  744.   }
  745.  
  746.   /* Pass 2: process columns. */
  747.   /* Note that we must descale the results by a factor of 8 == 2**3, */
  748.   /* and also undo the PASS1_BITS scaling. */
  749.  
  750.   dataptr = data;
  751.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  752.     /* Columns of zeroes can be exploited in the same way as we did with rows.
  753.      * However, the row calculation has created many nonzero AC terms, so the
  754.      * simplification applies less often (typically 5% to 10% of the time).
  755.      * On machines with very fast multiplication, it's possible that the
  756.      * test takes more time than it's worth.  In that case this section
  757.      * may be commented out.
  758.      */
  759.  
  760.     d0 = dataptr[DCTSIZE*0];
  761.     d1 = dataptr[DCTSIZE*1];
  762.     d2 = dataptr[DCTSIZE*2];
  763.     d3 = dataptr[DCTSIZE*3];
  764.     d4 = dataptr[DCTSIZE*4];
  765.     d5 = dataptr[DCTSIZE*5];
  766.     d6 = dataptr[DCTSIZE*6];
  767.     d7 = dataptr[DCTSIZE*7];
  768.  
  769.     /* Even part: reverse the even part of the forward DCT. */
  770.     /* The rotator is sqrt(2)*c(-6). */
  771.     if (d6) {
  772.     if (d4) {
  773.         if (d2) {
  774.         if (d0) {
  775.             /* d0 != 0, d2 != 0, d4 != 0, d6 != 0 */
  776.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  777.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  778.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  779.  
  780.             tmp0 = (d0 + d4) << CONST_BITS;
  781.             tmp1 = (d0 - d4) << CONST_BITS;
  782.  
  783.             tmp10 = tmp0 + tmp3;
  784.             tmp13 = tmp0 - tmp3;
  785.             tmp11 = tmp1 + tmp2;
  786.             tmp12 = tmp1 - tmp2;
  787.         } else {
  788.             /* d0 == 0, d2 != 0, d4 != 0, d6 != 0 */
  789.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  790.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  791.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  792.  
  793.             tmp0 = d4 << CONST_BITS;
  794.  
  795.             tmp10 = tmp0 + tmp3;
  796.             tmp13 = tmp0 - tmp3;
  797.             tmp11 = tmp2 - tmp0;
  798.             tmp12 = -(tmp0 + tmp2);
  799.         }
  800.         } else {
  801.         if (d0) {
  802.             /* d0 != 0, d2 == 0, d4 != 0, d6 != 0 */
  803.             tmp2 = MULTIPLY(d6, - FIX(1.306562965));
  804.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  805.  
  806.             tmp0 = (d0 + d4) << CONST_BITS;
  807.             tmp1 = (d0 - d4) << CONST_BITS;
  808.  
  809.             tmp10 = tmp0 + tmp3;
  810.             tmp13 = tmp0 - tmp3;
  811.             tmp11 = tmp1 + tmp2;
  812.             tmp12 = tmp1 - tmp2;
  813.         } else {
  814.             /* d0 == 0, d2 == 0, d4 != 0, d6 != 0 */
  815.             tmp2 = MULTIPLY(d6, -FIX(1.306562965));
  816.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  817.  
  818.             tmp0 = d4 << CONST_BITS;
  819.  
  820.             tmp10 = tmp0 + tmp3;
  821.             tmp13 = tmp0 - tmp3;
  822.             tmp11 = tmp2 - tmp0;
  823.             tmp12 = -(tmp0 + tmp2);
  824.         }
  825.         }
  826.     } else {
  827.         if (d2) {
  828.         if (d0) {
  829.             /* d0 != 0, d2 != 0, d4 == 0, d6 != 0 */
  830.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  831.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  832.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  833.  
  834.             tmp0 = d0 << CONST_BITS;
  835.  
  836.             tmp10 = tmp0 + tmp3;
  837.             tmp13 = tmp0 - tmp3;
  838.             tmp11 = tmp0 + tmp2;
  839.             tmp12 = tmp0 - tmp2;
  840.         } else {
  841.             /* d0 == 0, d2 != 0, d4 == 0, d6 != 0 */
  842.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  843.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  844.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  845.  
  846.             tmp10 = tmp3;
  847.             tmp13 = -tmp3;
  848.             tmp11 = tmp2;
  849.             tmp12 = -tmp2;
  850.         }
  851.         } else {
  852.         if (d0) {
  853.             /* d0 != 0, d2 == 0, d4 == 0, d6 != 0 */
  854.             tmp2 = MULTIPLY(d6, - FIX(1.306562965));
  855.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  856.  
  857.             tmp0 = d0 << CONST_BITS;
  858.  
  859.             tmp10 = tmp0 + tmp3;
  860.             tmp13 = tmp0 - tmp3;
  861.             tmp11 = tmp0 + tmp2;
  862.             tmp12 = tmp0 - tmp2;
  863.         } else {
  864.             /* d0 == 0, d2 == 0, d4 == 0, d6 != 0 */
  865.             tmp2 = MULTIPLY(d6, - FIX(1.306562965));
  866.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  867.  
  868.             tmp10 = tmp3;
  869.             tmp13 = -tmp3;
  870.             tmp11 = tmp2;
  871.             tmp12 = -tmp2;
  872.         }
  873.         }
  874.     }
  875.     } else {
  876.     if (d4) {
  877.         if (d2) {
  878.         if (d0) {
  879.             /* d0 != 0, d2 != 0, d4 != 0, d6 == 0 */
  880.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  881.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  882.  
  883.             tmp0 = (d0 + d4) << CONST_BITS;
  884.             tmp1 = (d0 - d4) << CONST_BITS;
  885.  
  886.             tmp10 = tmp0 + tmp3;
  887.             tmp13 = tmp0 - tmp3;
  888.             tmp11 = tmp1 + tmp2;
  889.             tmp12 = tmp1 - tmp2;
  890.         } else {
  891.             /* d0 == 0, d2 != 0, d4 != 0, d6 == 0 */
  892.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  893.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  894.  
  895.             tmp0 = d4 << CONST_BITS;
  896.  
  897.             tmp10 = tmp0 + tmp3;
  898.             tmp13 = tmp0 - tmp3;
  899.             tmp11 = tmp2 - tmp0;
  900.             tmp12 = -(tmp0 + tmp2);
  901.         }
  902.         } else {
  903.         if (d0) {
  904.             /* d0 != 0, d2 == 0, d4 != 0, d6 == 0 */
  905.             tmp10 = tmp13 = (d0 + d4) << CONST_BITS;
  906.             tmp11 = tmp12 = (d0 - d4) << CONST_BITS;
  907.         } else {
  908.             /* d0 == 0, d2 == 0, d4 != 0, d6 == 0 */
  909.             tmp10 = tmp13 = d4 << CONST_BITS;
  910.             tmp11 = tmp12 = -tmp10;
  911.         }
  912.         }
  913.     } else {
  914.         if (d2) {
  915.         if (d0) {
  916.             /* d0 != 0, d2 != 0, d4 == 0, d6 == 0 */
  917.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  918.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  919.  
  920.             tmp0 = d0 << CONST_BITS;
  921.  
  922.             tmp10 = tmp0 + tmp3;
  923.             tmp13 = tmp0 - tmp3;
  924.             tmp11 = tmp0 + tmp2;
  925.             tmp12 = tmp0 - tmp2;
  926.         } else {
  927.             /* d0 == 0, d2 != 0, d4 == 0, d6 == 0 */
  928.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  929.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  930.  
  931.             tmp10 = tmp3;
  932.             tmp13 = -tmp3;
  933.             tmp11 = tmp2;
  934.             tmp12 = -tmp2;
  935.         }
  936.         } else {
  937.         if (d0) {
  938.             /* d0 != 0, d2 == 0, d4 == 0, d6 == 0 */
  939.             tmp10 = tmp13 = tmp11 = tmp12 = d0 << CONST_BITS;
  940.         } else {
  941.             /* d0 == 0, d2 == 0, d4 == 0, d6 == 0 */
  942.             tmp10 = tmp13 = tmp11 = tmp12 = 0;
  943.         }
  944.         }
  945.     }
  946.     }
  947.  
  948.     /* Odd part per figure 8; the matrix is unitary and hence its
  949.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  950.      */
  951.     if (d7) {
  952.     if (d5) {
  953.         if (d3) {
  954.         if (d1) {
  955.             /* d1 != 0, d3 != 0, d5 != 0, d7 != 0 */
  956.             z1 = d7 + d1;
  957.             z2 = d5 + d3;
  958.             z3 = d7 + d3;
  959.             z4 = d5 + d1;
  960.             z5 = MULTIPLY(z3 + z4, FIX(1.175875602));
  961.             
  962.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  963.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  964.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  965.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  966.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  967.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  968.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  969.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  970.             
  971.             z3 += z5;
  972.             z4 += z5;
  973.             
  974.             tmp0 += z1 + z3;
  975.             tmp1 += z2 + z4;
  976.             tmp2 += z2 + z3;
  977.             tmp3 += z1 + z4;
  978.         } else {
  979.             /* d1 == 0, d3 != 0, d5 != 0, d7 != 0 */
  980.             z1 = d7;
  981.             z2 = d5 + d3;
  982.             z3 = d7 + d3;
  983.             z5 = MULTIPLY(z3 + d5, FIX(1.175875602));
  984.             
  985.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  986.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  987.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  988.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  989.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  990.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  991.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  992.             
  993.             z3 += z5;
  994.             z4 += z5;
  995.             
  996.             tmp0 += z1 + z3;
  997.             tmp1 += z2 + z4;
  998.             tmp2 += z2 + z3;
  999.             tmp3 = z1 + z4;
  1000.         }
  1001.         } else {
  1002.         if (d1) {
  1003.             /* d1 != 0, d3 == 0, d5 != 0, d7 != 0 */
  1004.             z1 = d7 + d1;
  1005.             z2 = d5;
  1006.             z3 = d7;
  1007.             z4 = d5 + d1;
  1008.             z5 = MULTIPLY(z3 + z4, FIX(1.175875602));
  1009.             
  1010.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  1011.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  1012.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  1013.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  1014.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  1015.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  1016.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  1017.             
  1018.             z3 += z5;
  1019.             z4 += z5;
  1020.             
  1021.             tmp0 += z1 + z3;
  1022.             tmp1 += z2 + z4;
  1023.             tmp2 = z2 + z3;
  1024.             tmp3 += z1 + z4;
  1025.         } else {
  1026.             /* d1 == 0, d3 == 0, d5 != 0, d7 != 0 */
  1027.             tmp0 = MULTIPLY(d7, - FIX(0.601344887)); 
  1028.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  1029.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  1030.             tmp1 = MULTIPLY(d5, - FIX(0.509795578));
  1031.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  1032.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  1033.             z5 = MULTIPLY(d5 + d7, FIX(1.175875602));
  1034.             
  1035.             z3 += z5;
  1036.             z4 += z5;
  1037.             
  1038.             tmp0 += z3;
  1039.             tmp1 += z4;
  1040.             tmp2 = z2 + z3;
  1041.             tmp3 = z1 + z4;
  1042.         }
  1043.         }
  1044.     } else {
  1045.         if (d3) {
  1046.         if (d1) {
  1047.             /* d1 != 0, d3 != 0, d5 == 0, d7 != 0 */
  1048.             z1 = d7 + d1;
  1049.             z3 = d7 + d3;
  1050.             z5 = MULTIPLY(z3 + d1, FIX(1.175875602));
  1051.             
  1052.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  1053.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  1054.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  1055.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  1056.             z2 = MULTIPLY(d3, - FIX(2.562915447));
  1057.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  1058.             z4 = MULTIPLY(d1, - FIX(0.390180644));
  1059.             
  1060.             z3 += z5;
  1061.             z4 += z5;
  1062.             
  1063.             tmp0 += z1 + z3;
  1064.             tmp1 = z2 + z4;
  1065.             tmp2 += z2 + z3;
  1066.             tmp3 += z1 + z4;
  1067.         } else {
  1068.             /* d1 == 0, d3 != 0, d5 == 0, d7 != 0 */
  1069.             z3 = d7 + d3;
  1070.             
  1071.             tmp0 = MULTIPLY(d7, - FIX(0.601344887)); 
  1072.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  1073.             tmp2 = MULTIPLY(d3, FIX(0.509795579));
  1074.             z2 = MULTIPLY(d3, - FIX(2.562915447));
  1075.             z5 = MULTIPLY(z3, FIX(1.175875602));
  1076.             z3 = MULTIPLY(z3, - FIX(0.785694958));
  1077.             
  1078.             tmp0 += z3;
  1079.             tmp1 = z2 + z5;
  1080.             tmp2 += z3;
  1081.             tmp3 = z1 + z5;
  1082.         }
  1083.         } else {
  1084.         if (d1) {
  1085.             /* d1 != 0, d3 == 0, d5 == 0, d7 != 0 */
  1086.             z1 = d7 + d1;
  1087.             z5 = MULTIPLY(z1, FIX(1.175875602));
  1088.  
  1089.             z1 = MULTIPLY(z1, FIX(0.275899379));
  1090.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  1091.             tmp0 = MULTIPLY(d7, - FIX(1.662939224)); 
  1092.             z4 = MULTIPLY(d1, - FIX(0.390180644));
  1093.             tmp3 = MULTIPLY(d1, FIX(1.111140466));
  1094.  
  1095.             tmp0 += z1;
  1096.             tmp1 = z4 + z5;
  1097.             tmp2 = z3 + z5;
  1098.             tmp3 += z1;
  1099.         } else {
  1100.             /* d1 == 0, d3 == 0, d5 == 0, d7 != 0 */
  1101.             tmp0 = MULTIPLY(d7, - FIX(1.387039845));
  1102.             tmp1 = MULTIPLY(d7, FIX(1.175875602));
  1103.             tmp2 = MULTIPLY(d7, - FIX(0.785694958));
  1104.             tmp3 = MULTIPLY(d7, FIX(0.275899379));
  1105.         }
  1106.         }
  1107.     }
  1108.     } else {
  1109.     if (d5) {
  1110.         if (d3) {
  1111.         if (d1) {
  1112.             /* d1 != 0, d3 != 0, d5 != 0, d7 == 0 */
  1113.             z2 = d5 + d3;
  1114.             z4 = d5 + d1;
  1115.             z5 = MULTIPLY(d3 + z4, FIX(1.175875602));
  1116.             
  1117.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  1118.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  1119.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  1120.             z1 = MULTIPLY(d1, - FIX(0.899976223));
  1121.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  1122.             z3 = MULTIPLY(d3, - FIX(1.961570560));
  1123.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  1124.             
  1125.             z3 += z5;
  1126.             z4 += z5;
  1127.             
  1128.             tmp0 = z1 + z3;
  1129.             tmp1 += z2 + z4;
  1130.             tmp2 += z2 + z3;
  1131.             tmp3 += z1 + z4;
  1132.         } else {
  1133.             /* d1 == 0, d3 != 0, d5 != 0, d7 == 0 */
  1134.             z2 = d5 + d3;
  1135.             
  1136.             z5 = MULTIPLY(z2, FIX(1.175875602));
  1137.             tmp1 = MULTIPLY(d5, FIX(1.662939225));
  1138.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  1139.             z2 = MULTIPLY(z2, - FIX(1.387039845));
  1140.             tmp2 = MULTIPLY(d3, FIX(1.111140466));
  1141.             z3 = MULTIPLY(d3, - FIX(1.961570560));
  1142.             
  1143.             tmp0 = z3 + z5;
  1144.             tmp1 += z2;
  1145.             tmp2 += z2;
  1146.             tmp3 = z4 + z5;
  1147.         }
  1148.         } else {
  1149.         if (d1) {
  1150.             /* d1 != 0, d3 == 0, d5 != 0, d7 == 0 */
  1151.             z4 = d5 + d1;
  1152.             
  1153.             z5 = MULTIPLY(z4, FIX(1.175875602));
  1154.             z1 = MULTIPLY(d1, - FIX(0.899976223));
  1155.             tmp3 = MULTIPLY(d1, FIX(0.601344887));
  1156.             tmp1 = MULTIPLY(d5, - FIX(0.509795578));
  1157.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  1158.             z4 = MULTIPLY(z4, FIX(0.785694958));
  1159.             
  1160.             tmp0 = z1 + z5;
  1161.             tmp1 += z4;
  1162.             tmp2 = z2 + z5;
  1163.             tmp3 += z4;
  1164.         } else {
  1165.             /* d1 == 0, d3 == 0, d5 != 0, d7 == 0 */
  1166.             tmp0 = MULTIPLY(d5, FIX(1.175875602));
  1167.             tmp1 = MULTIPLY(d5, FIX(0.275899380));
  1168.             tmp2 = MULTIPLY(d5, - FIX(1.387039845));
  1169.             tmp3 = MULTIPLY(d5, FIX(0.785694958));
  1170.         }
  1171.         }
  1172.     } else {
  1173.         if (d3) {
  1174.         if (d1) {
  1175.             /* d1 != 0, d3 != 0, d5 == 0, d7 == 0 */
  1176.             z5 = d1 + d3;
  1177.             tmp3 = MULTIPLY(d1, FIX(0.211164243));
  1178.             tmp2 = MULTIPLY(d3, - FIX(1.451774981));
  1179.             z1 = MULTIPLY(d1, FIX(1.061594337));
  1180.             z2 = MULTIPLY(d3, - FIX(2.172734803));
  1181.             z4 = MULTIPLY(z5, FIX(0.785694958));
  1182.             z5 = MULTIPLY(z5, FIX(1.175875602));
  1183.             
  1184.             tmp0 = z1 - z4;
  1185.             tmp1 = z2 + z4;
  1186.             tmp2 += z5;
  1187.             tmp3 += z5;
  1188.         } else {
  1189.             /* d1 == 0, d3 != 0, d5 == 0, d7 == 0 */
  1190.             tmp0 = MULTIPLY(d3, - FIX(0.785694958));
  1191.             tmp1 = MULTIPLY(d3, - FIX(1.387039845));
  1192.             tmp2 = MULTIPLY(d3, - FIX(0.275899379));
  1193.             tmp3 = MULTIPLY(d3, FIX(1.175875602));
  1194.         }
  1195.         } else {
  1196.         if (d1) {
  1197.             /* d1 != 0, d3 == 0, d5 == 0, d7 == 0 */
  1198.             tmp0 = MULTIPLY(d1, FIX(0.275899379));
  1199.             tmp1 = MULTIPLY(d1, FIX(0.785694958));
  1200.             tmp2 = MULTIPLY(d1, FIX(1.175875602));
  1201.             tmp3 = MULTIPLY(d1, FIX(1.387039845));
  1202.         } else {
  1203.             /* d1 == 0, d3 == 0, d5 == 0, d7 == 0 */
  1204.             tmp0 = tmp1 = tmp2 = tmp3 = 0;
  1205.         }
  1206.         }
  1207.     }
  1208.     }
  1209.  
  1210.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  1211.  
  1212.     dataptr[DCTSIZE*0] = (DCTELEM) DESCALE(tmp10 + tmp3,
  1213.                        CONST_BITS+PASS1_BITS+3);
  1214.     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp10 - tmp3,
  1215.                        CONST_BITS+PASS1_BITS+3);
  1216.     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp11 + tmp2,
  1217.                        CONST_BITS+PASS1_BITS+3);
  1218.     dataptr[DCTSIZE*6] = (DCTELEM) DESCALE(tmp11 - tmp2,
  1219.                        CONST_BITS+PASS1_BITS+3);
  1220.     dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(tmp12 + tmp1,
  1221.                        CONST_BITS+PASS1_BITS+3);
  1222.     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp12 - tmp1,
  1223.                        CONST_BITS+PASS1_BITS+3);
  1224.     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp13 + tmp0,
  1225.                        CONST_BITS+PASS1_BITS+3);
  1226.     dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(tmp13 - tmp0,
  1227.                        CONST_BITS+PASS1_BITS+3);
  1228.     
  1229.     dataptr++;            /* advance pointer to next column */
  1230.   }
  1231. }
  1232.  
  1233. #else
  1234.  
  1235.  
  1236. void
  1237. j_rev_dct_sparse (data, pos) 
  1238.      DCTBLOCK data;
  1239.      int pos;
  1240. {
  1241.   j_rev_dct(data);
  1242. }
  1243.  
  1244. void
  1245. j_rev_dct (data)
  1246.   DCTBLOCK data;
  1247. {
  1248.   INT32 tmp0, tmp1, tmp2, tmp3;
  1249.   INT32 tmp10, tmp11, tmp12, tmp13;
  1250.   INT32 z1, z2, z3, z4, z5;
  1251.   register DCTELEM *dataptr;
  1252.   int rowctr;
  1253.   SHIFT_TEMPS
  1254.  
  1255.   /* Pass 1: process rows. */
  1256.   /* Note results are scaled up by sqrt(8) compared to a true IDCT; */
  1257.   /* furthermore, we scale the results by 2**PASS1_BITS. */
  1258.  
  1259.   dataptr = data;
  1260.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  1261.     /* Due to quantization, we will usually find that many of the input
  1262.      * coefficients are zero, especially the AC terms.  We can exploit this
  1263.      * by short-circuiting the IDCT calculation for any row in which all
  1264.      * the AC terms are zero.  In that case each output is equal to the
  1265.      * DC coefficient (with scale factor as needed).
  1266.      * With typical images and quantization tables, half or more of the
  1267.      * row DCT calculations can be simplified this way.
  1268.      */
  1269.  
  1270.     if ((dataptr[1] | dataptr[2] | dataptr[3] | dataptr[4] |
  1271.      dataptr[5] | dataptr[6] | dataptr[7]) == 0) {
  1272.       /* AC terms all zero */
  1273.       DCTELEM dcval = (DCTELEM) (dataptr[0] << PASS1_BITS);
  1274.       
  1275.       dataptr[0] = dcval;
  1276.       dataptr[1] = dcval;
  1277.       dataptr[2] = dcval;
  1278.       dataptr[3] = dcval;
  1279.       dataptr[4] = dcval;
  1280.       dataptr[5] = dcval;
  1281.       dataptr[6] = dcval;
  1282.       dataptr[7] = dcval;
  1283.       
  1284.       dataptr += DCTSIZE;    /* advance pointer to next row */
  1285.       continue;
  1286.     }
  1287.  
  1288.     /* Even part: reverse the even part of the forward DCT. */
  1289.     /* The rotator is sqrt(2)*c(-6). */
  1290.  
  1291.     z2 = (INT32) dataptr[2];
  1292.     z3 = (INT32) dataptr[6];
  1293.  
  1294.     z1 = MULTIPLY(z2 + z3, FIX(0.541196100));
  1295.     tmp2 = z1 + MULTIPLY(z3, - FIX(1.847759065));
  1296.     tmp3 = z1 + MULTIPLY(z2, FIX(0.765366865));
  1297.  
  1298.     tmp0 = ((INT32) dataptr[0] + (INT32) dataptr[4]) << CONST_BITS;
  1299.     tmp1 = ((INT32) dataptr[0] - (INT32) dataptr[4]) << CONST_BITS;
  1300.  
  1301.     tmp10 = tmp0 + tmp3;
  1302.     tmp13 = tmp0 - tmp3;
  1303.     tmp11 = tmp1 + tmp2;
  1304.     tmp12 = tmp1 - tmp2;
  1305.     
  1306.     /* Odd part per figure 8; the matrix is unitary and hence its
  1307.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  1308.      */
  1309.  
  1310.     tmp0 = (INT32) dataptr[7];
  1311.     tmp1 = (INT32) dataptr[5];
  1312.     tmp2 = (INT32) dataptr[3];
  1313.     tmp3 = (INT32) dataptr[1];
  1314.  
  1315.     z1 = tmp0 + tmp3;
  1316.     z2 = tmp1 + tmp2;
  1317.     z3 = tmp0 + tmp2;
  1318.     z4 = tmp1 + tmp3;
  1319.     z5 = MULTIPLY(z3 + z4, FIX(1.175875602)); /* sqrt(2) * c3 */
  1320.     
  1321.     tmp0 = MULTIPLY(tmp0, FIX(0.298631336)); /* sqrt(2) * (-c1+c3+c5-c7) */
  1322.     tmp1 = MULTIPLY(tmp1, FIX(2.053119869)); /* sqrt(2) * ( c1+c3-c5+c7) */
  1323.     tmp2 = MULTIPLY(tmp2, FIX(3.072711026)); /* sqrt(2) * ( c1+c3+c5-c7) */
  1324.     tmp3 = MULTIPLY(tmp3, FIX(1.501321110)); /* sqrt(2) * ( c1+c3-c5-c7) */
  1325.     z1 = MULTIPLY(z1, - FIX(0.899976223)); /* sqrt(2) * (c7-c3) */
  1326.     z2 = MULTIPLY(z2, - FIX(2.562915447)); /* sqrt(2) * (-c1-c3) */
  1327.     z3 = MULTIPLY(z3, - FIX(1.961570560)); /* sqrt(2) * (-c3-c5) */
  1328.     z4 = MULTIPLY(z4, - FIX(0.390180644)); /* sqrt(2) * (c5-c3) */
  1329.     
  1330.     z3 += z5;
  1331.     z4 += z5;
  1332.     
  1333.     tmp0 += z1 + z3;
  1334.     tmp1 += z2 + z4;
  1335.     tmp2 += z2 + z3;
  1336.     tmp3 += z1 + z4;
  1337.  
  1338.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  1339.  
  1340.     dataptr[0] = (DCTELEM) DESCALE(tmp10 + tmp3, CONST_BITS-PASS1_BITS);
  1341.     dataptr[7] = (DCTELEM) DESCALE(tmp10 - tmp3, CONST_BITS-PASS1_BITS);
  1342.     dataptr[1] = (DCTELEM) DESCALE(tmp11 + tmp2, CONST_BITS-PASS1_BITS);
  1343.     dataptr[6] = (DCTELEM) DESCALE(tmp11 - tmp2, CONST_BITS-PASS1_BITS);
  1344.     dataptr[2] = (DCTELEM) DESCALE(tmp12 + tmp1, CONST_BITS-PASS1_BITS);
  1345.     dataptr[5] = (DCTELEM) DESCALE(tmp12 - tmp1, CONST_BITS-PASS1_BITS);
  1346.     dataptr[3] = (DCTELEM) DESCALE(tmp13 + tmp0, CONST_BITS-PASS1_BITS);
  1347.     dataptr[4] = (DCTELEM) DESCALE(tmp13 - tmp0, CONST_BITS-PASS1_BITS);
  1348.  
  1349.     dataptr += DCTSIZE;        /* advance pointer to next row */
  1350.   }
  1351.  
  1352.   /* Pass 2: process columns. */
  1353.   /* Note that we must descale the results by a factor of 8 == 2**3, */
  1354.   /* and also undo the PASS1_BITS scaling. */
  1355.  
  1356.   dataptr = data;
  1357.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  1358.     /* Columns of zeroes can be exploited in the same way as we did with rows.
  1359.      * However, the row calculation has created many nonzero AC terms, so the
  1360.      * simplification applies less often (typically 5% to 10% of the time).
  1361.      * On machines with very fast multiplication, it's possible that the
  1362.      * test takes more time than it's worth.  In that case this section
  1363.      * may be commented out.
  1364.      */
  1365.  
  1366. #ifndef NO_ZERO_COLUMN_TEST
  1367.     if ((dataptr[DCTSIZE*1] | dataptr[DCTSIZE*2] | dataptr[DCTSIZE*3] |
  1368.      dataptr[DCTSIZE*4] | dataptr[DCTSIZE*5] | dataptr[DCTSIZE*6] |
  1369.      dataptr[DCTSIZE*7]) == 0) {
  1370.       /* AC terms all zero */
  1371.       DCTELEM dcval = (DCTELEM) DESCALE((INT32) dataptr[0], PASS1_BITS+3);
  1372.       
  1373.       dataptr[DCTSIZE*0] = dcval;
  1374.       dataptr[DCTSIZE*1] = dcval;
  1375.       dataptr[DCTSIZE*2] = dcval;
  1376.       dataptr[DCTSIZE*3] = dcval;
  1377.       dataptr[DCTSIZE*4] = dcval;
  1378.       dataptr[DCTSIZE*5] = dcval;
  1379.       dataptr[DCTSIZE*6] = dcval;
  1380.       dataptr[DCTSIZE*7] = dcval;
  1381.       
  1382.       dataptr++;        /* advance pointer to next column */
  1383.       continue;
  1384.     }
  1385. #endif
  1386.  
  1387.     /* Even part: reverse the even part of the forward DCT. */
  1388.     /* The rotator is sqrt(2)*c(-6). */
  1389.  
  1390.     z2 = (INT32) dataptr[DCTSIZE*2];
  1391.     z3 = (INT32) dataptr[DCTSIZE*6];
  1392.  
  1393.     z1 = MULTIPLY(z2 + z3, FIX(0.541196100));
  1394.     tmp2 = z1 + MULTIPLY(z3, - FIX(1.847759065));
  1395.     tmp3 = z1 + MULTIPLY(z2, FIX(0.765366865));
  1396.  
  1397.     tmp0 = ((INT32) dataptr[DCTSIZE*0] + (INT32) dataptr[DCTSIZE*4]) << CONST_BITS;
  1398.     tmp1 = ((INT32) dataptr[DCTSIZE*0] - (INT32) dataptr[DCTSIZE*4]) << CONST_BITS;
  1399.  
  1400.     tmp10 = tmp0 + tmp3;
  1401.     tmp13 = tmp0 - tmp3;
  1402.     tmp11 = tmp1 + tmp2;
  1403.     tmp12 = tmp1 - tmp2;
  1404.     
  1405.     /* Odd part per figure 8; the matrix is unitary and hence its
  1406.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  1407.      */
  1408.  
  1409.     tmp0 = (INT32) dataptr[DCTSIZE*7];
  1410.     tmp1 = (INT32) dataptr[DCTSIZE*5];
  1411.     tmp2 = (INT32) dataptr[DCTSIZE*3];
  1412.     tmp3 = (INT32) dataptr[DCTSIZE*1];
  1413.  
  1414.     z1 = tmp0 + tmp3;
  1415.     z2 = tmp1 + tmp2;
  1416.     z3 = tmp0 + tmp2;
  1417.     z4 = tmp1 + tmp3;
  1418.     z5 = MULTIPLY(z3 + z4, FIX(1.175875602)); /* sqrt(2) * c3 */
  1419.     
  1420.     tmp0 = MULTIPLY(tmp0, FIX(0.298631336)); /* sqrt(2) * (-c1+c3+c5-c7) */
  1421.     tmp1 = MULTIPLY(tmp1, FIX(2.053119869)); /* sqrt(2) * ( c1+c3-c5+c7) */
  1422.     tmp2 = MULTIPLY(tmp2, FIX(3.072711026)); /* sqrt(2) * ( c1+c3+c5-c7) */
  1423.     tmp3 = MULTIPLY(tmp3, FIX(1.501321110)); /* sqrt(2) * ( c1+c3-c5-c7) */
  1424.     z1 = MULTIPLY(z1, - FIX(0.899976223)); /* sqrt(2) * (c7-c3) */
  1425.     z2 = MULTIPLY(z2, - FIX(2.562915447)); /* sqrt(2) * (-c1-c3) */
  1426.     z3 = MULTIPLY(z3, - FIX(1.961570560)); /* sqrt(2) * (-c3-c5) */
  1427.     z4 = MULTIPLY(z4, - FIX(0.390180644)); /* sqrt(2) * (c5-c3) */
  1428.     
  1429.     z3 += z5;
  1430.     z4 += z5;
  1431.     
  1432.     tmp0 += z1 + z3;
  1433.     tmp1 += z2 + z4;
  1434.     tmp2 += z2 + z3;
  1435.     tmp3 += z1 + z4;
  1436.  
  1437.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  1438.  
  1439.     dataptr[DCTSIZE*0] = (DCTELEM) DESCALE(tmp10 + tmp3,
  1440.                        CONST_BITS+PASS1_BITS+3);
  1441.     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp10 - tmp3,
  1442.                        CONST_BITS+PASS1_BITS+3);
  1443.     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp11 + tmp2,
  1444.                        CONST_BITS+PASS1_BITS+3);
  1445.     dataptr[DCTSIZE*6] = (DCTELEM) DESCALE(tmp11 - tmp2,
  1446.                        CONST_BITS+PASS1_BITS+3);
  1447.     dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(tmp12 + tmp1,
  1448.                        CONST_BITS+PASS1_BITS+3);
  1449.     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp12 - tmp1,
  1450.                        CONST_BITS+PASS1_BITS+3);
  1451.     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp13 + tmp0,
  1452.                        CONST_BITS+PASS1_BITS+3);
  1453.     dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(tmp13 - tmp0,
  1454.                        CONST_BITS+PASS1_BITS+3);
  1455.     
  1456.     dataptr++;            /* advance pointer to next column */
  1457.   }
  1458. }
  1459.  
  1460.  
  1461. #endif
  1462.